library(plyr)
library(robust)

f = oda ~ p_polity2 + pop
f1 = update(f, '.~.+alliance*independence')
f2 = update(f, '.~.+trade*independence')
f3 = update(f, '.~.+infant_mortality*independence')

cols = c('agency', 'iso3c', 'year', all.vars(f), 'trade', 'alliance', 'infant_mortality', 'independence')
z = na.omit(dat[,cols])

f = function(i) colMeans(i[,4:10])
k = ddply(z, c('agency', 'iso3c'), f)

mod1 = lm(f1, k)
mod2 = lm(f2, k)
mod3 = lm(f3, k)

mod4 = lmRob(f1, data=k, control=lmRob.control(mxr = 500, mxf = 500, mxs = 500, seed=1024))
mod5 = lmRob(f2, data=k, control=lmRob.control(mxr = 500, mxf = 500, mxs = 500, seed=1024))
mod6 = lmRob(f3, data=k, control=lmRob.control(mxr = 500, mxf = 500, mxs = 500, seed=1024))

se = lapply(list(mod1, mod2, mod3), sandwich)
se = lapply(se, function(i) sqrt(diag(i)))
